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Transient Water Age Distributions in Environmental Flow 
Systems: The Time-Marching Laplace Transform Solution 
Technique. 

F. J. Cornaton^ 

Abstract. Environmental fluid circulations are very often characterized by analyzing 

the fate and behavior of natural and anthropogenic tracers. Among these tracers, age 
is taken as an ideal tracer which can yield interesting diagnoses, as for example the char- 
acterization of the mixing and renewal of water masses, of the fate and mixing of con- 
taminants, or the calibration of hydro- dispersive parameters used by numerical models. 
Such diagnoses are of great interest in atmospheric and ocean circulation sciences, as well 
in surface and subsurface hydrology. The temporal evolution of groundwater age and its 
frequency distributions can display important changes as flow regimes vary due to nat- 
ural change in climate and hydrologic conditions and/or human induced pressures on the 
resource to satisfy the water demand. Groundwater age being nowadays frequently used 
to investigate reservoir properties and recharge conditions, special attention needs to be 
put on the way this property is characterized, would it be using isotopic methods or math- 
ematical modelling. Steady-state age frequency distributions can be modelled using stan- 
dard numerical techniques, since the general balance equation describing age transport 
under steady-state flow conditions is exactly equivalent to a standard advection-dispersion 
equation. The time-dependent problem is however described by an extended transport 
operator that incorporates an additional coordinate for water age. The consequence is 
that numerical solutions can hardly be achieved, especially for real 3-D applications over 
large time periods of interest. A novel algorithm for solving the age distribution prob- 
lem under time-varying flow regimes is prc^scmted and, for some specific configurations, 
extended to the problem of generalized component exposure time. The algorithm com- 
bines the Laplace Transform technique applied to the age (or exposure time) coordinate 
with standard time-marching schemes. The method is validated and illustrated using an- 
alytical and numerical solutions considering 1-D, 2-D and 3-D theoretical groundwater 
flow domains. 



1. Introduction 

A comprehensive analysis of flow processes and advcctivc- 
dispersive-difi^usive transport of dissolved constituents in 
natural fluid circulations is an important challenge in Earth 
and environmental sciences. An approach that has become 
very popular consists in using water age as an ideal tracer to 
tag fluid masses and estimate associated timescales {Deleer- 
snijder et al. [2010]). Such timescales can lead to very help- 
ful diagnoses that apply to interdisciplinary environmental 
studies, like in atmospheric and ocean circulation sciences 
or surface and subsurface hydrology. 

The characterization of water age characteristics in environ- 
mental flow systems is linked to the estimation of recharge 
patterns and variations, fluid flow dynamics over possibly 
various time scales, advcctive-dispcrsivc/diffusivc transport 
in heterogeneous circulations, and to some extend isotope 
geochemistry, fractionation, and interfacial reactive chem- 
istry {Glynn and Plummer [2005]; Ginn et al. [2009]). In 
particular the characterization of groundwater age distribu- 
tions is of prime interest as it can act as an indicator of 
aquifer contamination and vulnerability, explain the modes 
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of recharge of aquiferous reservoirs, or identify the diffu- 
sive exchanges between mobile and immobile flow regions. 
The literature offers a wide range of studies focusing on 
various aspects of water age and its conceptual representa/- 
tions. However because different mathematical and concej>- 
tual foundations are used, a standardization of the results 
is still diflicult to be achieved. 

Age is commonly inferred from concentrations of single iso- 
topes/tracers (for exhaustive reviews of groundwater age 
dating methods e.g. see Clarke and Fritz [2007]; Kazemi et 
al. [2006]), or from multiple isotope/tracer concentrations 
that together represent different components of age (e.g. 
sec Gorcho Alvarado et al. [2007]; Troldborg et al. [2008]; 
Larocque et al. [2009]; Lavastre et al. [2010]). This concept 
is however very misleading because concentration-based ages 
arc simply apparent ages that do not a priori relate to the 
mean age of a water sample {Sanford [2011]). Moreover, the 
use of environmental and anthropogenic tracers to make in- 
ferences on water age goes together with the use of lumped- 
parameter models that assume specific age distributions re- 
gardless the complexity of the flow domain. 
A relatively recent representation of age and its distri- 
butions is the one that relies on the use of distributed 
flow models. Accounting for the relevant transport char- 
acteristics affecting age distributions however requires the 
use of adequate governing transport equations. Regard- 
ing this point, new advances have permitted to consolidate 
the unification of the foundations on age fate by means of 
a physically-based derivation of a transport equation for 
mean age (Goode [1996]), percentile of age and moments 
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of age {Varm and Carrera [1998]), and later on, after the 
early works of Campana [1987], of equations describing the 
complete age frequency distribution in the context of ocean 
circulations modelling {Deleersmjder et al. [2001]; Delhez 
and Deleersnijder [2002]) and subsurface hydrology [Ginn 
[1999, 2000a, b]; Cornaton [2003]; Cornaton and Perrochet 
[2006a, b, 2007]; Gmn [2007]; Woolfenden and Gmn [2009]). 
Putting aside the way age is inferred (or modelled), the tran- 
sient nature of flow systems and its implication on the tran- 
sient nature of age distributions have never been successfully 
addressed, particularly for large space- and time-scales and 
real-site applications. This aspect is however of high impor- 
tance since the actual situations at which age properties are 
measured are a result of past-to-present perturbations {Ginn 
et al. [2009]), at short and large time-scales. Bethke and 
Johnson [2008] correctly pointed out that changing the field 
of groundwater age dating is a means for thinking about 
groundwater age in a new way. We can additionally ar- 
gue that the way age is modelled is also to be refined. In 
particular the time-dependency of age distributions play a 
fundamental role in the interpretation of age and tracer 
data. Transient groundwater age distributions are for exam- 
ple a result of natural transient hydrologic conditions and 
of human-induced modifications on the water cycles by ar- 
tificial withdraw and recharge, but they can also originate 
from more specific phenomena as the temporal dependency 
of fluid flow to fluid density and viscosity (e.g. in thermoha- 
line problems). When large time-scales are considered, the 
hydrodynamic conditions of a system are changing in time, 
in relation to the evolutions of geomorphology, system inter- 
nal structure, climate change and thus recharge conditions. 
Any modification induces changes in the age distributions. 
This article presents the outcomes of research works on tran- 
sient age distributions and their numerical solutions. We 
present a novel numerical formulation that is able to handle 
this problem in 3-D discretized domains of arbitrary com- 
plexity, and that applies to the most relevant numerical in- 
tegration techniques such as the finite element, finite vol- 
ume, and finite differences methods. The method is also 
particularly suited to fulfill the needs related to the mod- 
eling of hydrologic evolutions over very large time periods. 
In the first section we summarize the mathematical founda- 
tions of transient age and constituent exposure time density 
equations. In the second section we make a short review of 
existing numerical integration schemes and present a novel 
numerical formulation. The characteristics of the proposed 
algorithms are analyzed and discussed, and eventually vali- 
dated and illustrated in the last section using analytical and 
numerical solutions. 

2. IVlathematical Formulations 

This section summarizes the theoretical basics for mod- 
eling water age and constituent exposure time distributions 
under transient conditions of fluid flow. 



be the physical coordinates for the case n — 3, and t be the 
clock-time, a point (a;, t) in the physical distance and time 
dimensions lies only in the 4-D space R" x , while a point 
{x,t,T) lies in the augmented 5-D space R""*"^ x IS}. 
The equation governing the distribution of constituent a ex- 
posure time density pa = pci{x, t, r) within a flow domain 
is given by Ginn [1999]: 



dt 



dVapa 

^ 5 = ^° 

OT 



(1) 



in which = Ja a + Ja d the total mass density flux that 
accounts for advective flux ^ = vp^ and diffusive and 
pore-scaled dispersive flux The physical velocity field 

v(a;, t) derives from a pre-solution of an appropriate fluid 
flow equation. We use the Fickian constitutive theory for 
diffusion and dispersion in order to express ^ as a po- 
tential in mass density, d ~ "DVp^,, with Yi[x,t) being 
a conventional dispersion-diffusion tensor. The total mass 
density flux is then expressed as = vp^ — DVpu,. The 
exposure-time velocity Va is the rate of change of a mate- 
rial point location on the exposure-time direction r. The 
reaction term represents the material mass density rate 
of transformation. 



Transient-state distribution p {x, t, r) 




fnitiai distribution 



PM, 0, r) 



Figure 1. Illustration of a transient distribution of 
water age as deduced from an oscillating uniform one- 
dimensional velocity fleld. 



2.1. Equations for the Distribution of Constituent 
Exposure Time and Water age 

The time spent by constituent molecules undergoing 
transport in a given flow domain (e.g. water, solutes, sus- 
pended colloids) can be defined as an independent dimension 
over which elemental masses are distributed, just as they are 
distributed over a spatial dimension (see Gmn [1999]; Del- 
hez and Deleersnijder [2002]). Following Ginn [1999], the 
exposure time is defined as this independent dimension, and 
it is convenient to define the space over which such density 
distributions evolve as the n + 1 dimensional space, that is 
elaborated by augmenting the n physical dimensions with 
one real independent dimension, say r, that represents an 
exposure time used by constituent molecules in a given reser- 
voir. This additional dimension is a continuous measure or- 
thogonal to the remaining coordinates. Letting x — {x, y, z) 



For the case of water age [a = w), the density pw ~ 
pw{x, t, t) is the space- and time-dependent bulk- volumetric 
mass density of water distributed over age r, and the expo- 
sure time velocity in the r-direction is simply 1. Conse- 
quently, the equation governing water age density simplifies 



dpu, 
dt 



dpu 
dr 



(2) 



With equation (2), age distributes by the effect of the ad- 
vection at a unit velocity of water mass density along the 
age axis, and the age variation is expressed by the term . 
Age mixing is classically ruled by the variations of fluid ve- 
locity V and by the processes of dispersive/diffusive mixing 
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represented by the D-tensor. A theoretical density distribu- 
tion of water age density distribution pw{x,t,T) at a given 
position X is displayed in Figure 1. The distribution under- 
goes a temporal evolution from an initial situation that may 
correspond to an equilibrated (or steady) state. 

2.2. The Case of Groundwater Age 

The specific context of subsurface flow implies fluid mo- 
tion through heterogeneous porous media of porosity (or 
mobile water content) 9{x,t). The flow equation yielding 
the distribution of fluid flux <i{x,t) — 9{x,t)v{x,t) can be 
written using the modified form of Richards' equation, which 
describes the three-dimensional transient groundwater flow 
in a variably-saturated porous matrix: 



— -h V • q = gi - go 



(3) 



in which the fluid flux vector is q = —krKX/H with 
H{x, t) = h{x, t) + z being the hydraulic head and h{x, t) 
the pressure head, and where K is the saturated tensor of 
hydraulic conductivity, kr = kr{Sw) is the relative perme- 
ability of the medium with respect to the degree of water 
saturation = with 9^ being the saturated water con- 
tent. The terms qi and go are fluid source and sink terms, 
respectively. 

Using the definition of the bulk-volumetric water mass 
density pyj{x,t,T) = 9{x,t)g{x,t,T), with g{x,t,T) being 
the frequency distribution of groundwater ages, and formu- 
lating the dispersive component of mass fiux as — 6DVp, 
equation (2) transforms into 



deg 

dt 



= -V • (qp - ^DVfl) 



dr 



+ r 



(4) 



in which the term r = qi5{r) — qog accounts for the diver- 
gence of fiuid fiux in equation (3). Equation (4) is subject 
to the following initial and boundary conditions: 



g{x,Q,T) = fo{x,T) in fi 



g{x,t,T) = 5{t) 
-9DVg-n= h{t,T 
J ■ n = q(5(r) • n 
J • n = 



on r_,i 
) on r_,2 
on r_,3 

on To 



(5a) 

(5b) 
(5c) 
(5d) 

(5e) 



where fo{x,T) is a given function representing the (initial) 
age distribution function at time t = 0, n is a normal out- 
ward unit vector, and where the total age flux is defined 
by J = qp — 9'D'Vg. The boundary portions Fq and r_ de- 
note the no-flow and inflowing boundaries of O, respectively, 
and the indices 1, 2 and 3 refer to boundary portions along 
which Diriclilct, Neumann and Cauchy age conditions are 
prescribed, respectively. The zero-age boundary condition is 
mathematically represented by the delta distribution 5(r), 
which ensures a pure unit and instantaneous pulse input 
along the recharge boundary portions r_. Since the age- 
coordinate is totally independent on the time-coordinate, 
this condition is simply assigned in a natural way. Even if 
both the Dirichlet and Cauchy conditions (5b) and (5d) can 
be combined, the Cauchy condition appears to be the one 
that is the most meaningful from a physical point of view, 
because it prevents water molecules from moving upstream 
and exiting the system by the inlet limits (it is homogeneous 
for T > 0). The Dirichlet- type condition is also homogeneous 
for r > 0; however it does not prevent the backward move- 
ment of water particles at inflowing limits when dispersion 
is significant. Moreover, the use of the Dirichlet-type con- 
dition can yield unphysical results when some parts of the 
inlet boundaries become inactive (i.e. during phases when 



recharge at the boundaries is not present). The Cauchy con- 
dition is active as long as inflow rates are non-zero (q • n 
in equation (5d)), and consequently it can be handled natu- 
rally without any special care, while the Dirichlet condition 
implies the controlled action of deactivating the condition 
where and when inflow limits stop producing, and of reac- 
tivating it as soon as these limits start to produce water 
again. The Neumann-type condition (5c) can also be used 
at particular boundary portions where the production of age 
mass is driven by diffusion only. The function /2 represents 
a purely diffusive production of age mass along r2 boundary 
portions. 

Particular cases for which the internal production of wai- 
ter is associated to more specific, not necessary equal to 
zero distributions of age r/ will require the volumetric fiuid 
source intensity qi to be associated to an age distribu- 
tion gi{x,t,T) — <5(r — Ti{t)), being itself potentially time- 
dependent, such that the age mass source will be distributed 
as qi{x,t)6{T - Ti{t)). 

The equation describing the evolution of mean age is ob- 
tained by taking the first age moment of equation (4), since 
the mean age a{x, t) of water molecules in a REV is ob- 
tained by taking the first age moment of the age density 
distribution function g{x,t,T): 



a{x, t)= [ 
Jo 



Tg{x, t, T)dT 



(6) 



Taking the first order age moment of equations (4) and (5) 
yields the following boundary value problem for mean age: 



d9a 



-V • (qa - 6'DVo) - goa + 6 



a{x, 0) = ao{x) 
a{x, t) = on 
-6>DVo- n = /2(t) 



in f2 

r-.i 

on r_,2 







r-.3 u r, 







(7a) 

(7b) 
(7c) 
(7d) 

(7e) 

where mean age flux is deflned as J a = qa—9D\/a. Eq. (7a) 
corresponds to the equation derived by {Goode [1996]), in 
which mean age is continuously generated during water flow 
since porosity is acting as a source term. This source term 
indicates that water is aging one unit per unit time, in av- 
erage. Equation (7a) is exactly equivalent to a standard 
advcction-dispcrsion equation, and can be solved by any 
standard code that can jiroducc advective-dispcrsivc solute 
transport solutions. The only i)articularity resides in a mass 
source term that needs to be distributed as porosity is. 
Higher-order moment equations can also be derived from 
equation (4) in analogy to the development yielding equa^ 
tion (7a) to eventually recover the moment equations derived 
by Vami and Carrera [1998] or Delhez and Deleersnijder 
[2002]. 

3. Numerical Solutions of the Transient 
Age Distribution Equation 

In this section we first make a short review of existing 
numerical schemes yielding numerical solution of the tran- 
sient age distribution problem, and discuss their limitations. 
A novel integration scheme is then presented and analyzed 
in terms of its adequacy for handling problems of arbitrary 
complexity and configuration. 

Steady-state models of flow regimes involve steady-state 
age distributions, with which important simplifications are 
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possible since the steady age distribution satisfies the time- 
translation invariance g{x,t,t — r) = g{x,0, —t), Vt £ 
[0, . . . , oo[ {Haine et al. [2008]). This property means that 
steady distributions only depend on the age coordinate t. 

As a consequence, the steady-state form of the age distri- 
bution equation (4) is exactly equivalent to an advcction- 
disporsion equation of transient nature, and can be solved 
using standard advcction-dispcrsion solvers (sec Cornaton 
[2003] , Cornaton and Perrochet [2006a] , Cornaton and Per- 
rochet [2006b]): 

^ = -V-(q<7-eDVfl) + r (8) 

The only technical difficulty for solving equation (8) is re- 
lated to the type of boundary condition required to prop- 
erly describe the zero-age condition. This zero-age condi- 
tion is assumed to be Dirac along inflow boundaries, and is 
classically modelled by means of a Dirichlct-typo condition 
g(x,T) — 5{t) and/or a Cauchy-typc mass flux condition 
J ■ n — q(x)^(r) ■ n. In order to properly handle these con- 
ditions, Cornaton [2003] proposed to make use of Laplace 
transforms (since the function S{t) transforms into 1) and of 
the LTG formalism {Sudicky [1989]) to obtain finite element 
solutions. 

Unsteady hydrodynamic conditions involve unsteady age 
distributions with which the time-translation invariance in- 
herent in equation (8) is not valid anymore, such that age 
is now distributed along two independent temporal dimen- 
sions: clock-time t and age r. This fact means that an age 
distribution for each clock-time and field point has to be 
considered and, as a consequence, the 5-D transient age dis- 
tribution equation presents severe technical difficulties in its 
numerical resolution in terms of computational burden and, 
above all, in terms of memory usage. In fact, the additional 
dimension induced by the age dimension requires the storage 
of a complete distribution (/(r) at any computational point 
of the system, and for all considered simulation times. 
To figure out the enormous computational requirements in- 
volved by the solution of the transient age distribution equa- 
tion, assume that we wish to represent g{x, t, r) as a discrete 
matrix that is discrctized with Nx points (or cells) over field 
points X, Nt points over clock-time t, and N-r points over 
age T. On the one hand, the storage requirements do not 
depend on the discretization method and are directly pro- 
portional to NxNtNr. On the other hand, the computer 
time required to compute the matrix does depend on the 
choice of the method. Modern numerical models can gener- 
ally handle up to iV^ = 0(10""^) and Nt ~ = 0(10^"^). 
The resulting matrix therefore has O(10^''~'^^) elements 
and is obviously impossible to store. 

Ginn [1999] proposed to treat the first-order term |^ in 
analogy to a physical advective process, thus grouping this 
term with the advection term V • qg'. That is, the first- 
order terms of equation (4) are taken together as V • q^g, 
where = {q, 1} is the velocity vector augmented by the 
age dimension, which presents the important property that 
its divergence in the (n + l)-dimensions space is given by 
the divergence in the standard Cartesian space of dimen- 
sion n {Ginn [1999]). A direct consequence of this for- 
mulation is that the dispersion tensor D also needs to be 
augmented with another dimension by adding a column 
and a row of zero entries (Do = {D,0}). Using the aug- 
mented velocity vector and dispersion tensor, equation (4) 
becomes equivalent to an advection-dispersion equation, 
= —V • {q^g — 6DaV(/) + r, which can be solved for 
n < 3 using conventional solute transport codes. The obvi- 
ous restriction of such a treatment of the age coordinate is 



that the physical dimension n cannot be bigger than 2, that 
is, at best a 3-D Cartesian domain can be used to solve for 
equation (4) in a 2-D spatial representation of a reservoir. 
Another difficulty is related to the need of a preUminary 
(optimal) definition for the size and the discretization of 
the physical mesh representing the age dimension. Since 
the variations in time of the ages are not known prior to 
calculation, a large enough grid has to be designed in or- 
der to capture the increase in age in time and a uniform 
discretization in the age direction for each physical mesh 
point needs to be defined (i.e. a uniform grid in the age 
direction). These aspects are quite restrictive and no ideal 
solution seems to exist. The technique is useful to provide 
validation examples for other approaches, but is obviously 
not practical for real-site applications. 

In the context of the modeling of oceanic circulations and 
surface water bodies, Delhez and Deleersnijder [2002] also 
proposed a numerical solution approach, but their scheme 
has never proved to be efficient to handle large domains and 
time-scales. As already discussed in the previous section, 
the term in equation (4) can be considered as an advec- 
tion term towards larger ages and, hence, discrctized with 
the techniques usually applied to advection terms. Equa- 
tion (4) can therefore be solved as an evolution equation in 
the 5-D space W^~^^ x M}, when using an operator splitting 
technique to handle separately the advection terms in the 
physical space and the age space {Delhez and Deleersnijder 
[2002]). When finite element, finite volume, or finite dif- 
ferences techniques are used, such a resolution involves not 
only the discretization of equation (4) on a usual 3-D mesh 
or grid but also in the age dimension by splitting the age do- 
main T € [0, Tmax] iD-to an appropriate number of age classes 
from minimum age to a maxirrmm age Tmax- Delhez and 
Deleersnijder [2002] propose to avoid the discretization of 
the ageing term and associated errors by introducing the 
change of variables (t, r) {t = t,T = t — r), which yields 
^ + ^ + ■■■ ^ +■■■' and thus removes the dif- 
ferential terms in the r— direction. The resolution of equa.- 
tion (4) finally corresponds to solving a number of classical 
uncoupled advection-dispersion equations corresponding to 
the different considered discrete age classes. As opposed to 
the mcsli-bascd method of Ginn [1999], the operator split- 
ting technique of Delhez and Deleersnijder [2002] has the 
advantage of being applicable to 3-D spaces, but it presents 
the disadvantage that a fixed number of age classes and a 
maximum allowed age have to be preliminary defined. The 
number of age classes can only be limited due to obvious 
memory storage limits, and this method may not operate 
with cases for which time-frames bigger than the maximum 
age present in the system have to be considered. 

3.1. The Time-Marching Laplace Transform Solution 
Technique (TMLT) 

In this section, a novel integration scheme is presented. 
The main aspect of the scheme is that a reduction in di- 
mension is achieved by mathematical transformation of the 
age dimension. The method combines the formalism of the 
LTG technique {Sudicky [1989]) and the techniques of time- 
marching classically used with mass transport solutions. 
The basic principle relies on the fact that age r and clock- 
time t are totally independent quantities, and that the rate 
of ageing does not depend on time t, since it is always one per 
unit time. In order to be relieved from the problem induced 
by the existence of an advection term in the age coordinate 
(see the two previous sections), and for the sake of reduc- 
ing the 5-dimensional system to a 4-dimensional system, a 
linear integral transform is applied to the age dimension r. 
The function g{x, t, r) being locally integrable on [0, oo), we 
can make use of an unilateral Laplace transform since the 
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necessary condition for its existence is fulfilled. 

When a Laplace transform is applied to the 5-D tran- 
sient age distribution equation (4) in the r— dimension, one 
obtains the following transformed transient age equation: 

^ = -V ■ (q.g - enVg) - sOg + mo + f (9) 

where g = g{x, t, s) = L{g{x, t, r), r, s} = e'"^ g{x, t, T)dT 
is the transformed state of the function g, with a £ 
denoting the complex Laplace variable and L{} the for- 
ward Laplace transformation operator, and where the trans- 
formed reaction term is f = qi{x,t) — qo{x,t)g. The zero- 
order source term mo = mo(a;, t) = 0{x, t)g{x, t, r = 0) orig- 
inates from the definition of the Laplace transform of the 
derivative of a function, i.e. L{^^^^|^^^, r, s} = sg{x, t, s) — 
g{x, t, 0). Equation (9) has the form of a standard transient- 
state advection-dispersion equation with first-order decay 
term (term —s9g). The complex function g{x,t,s) is thus 
transported in space and time. 

The full set of transformed initial and boundary condi- 
tions yielding solution of equation (9) can be formalized by: 

^(a;, 0, s) = /o(a;, s) in Q (10a) 

g{x,t,s) = l on Fi (10b) 

-O'DVg-n^ }2{t,s) on Fa (10c) 

J • n = on Fo (lOd) 

J ■ n = q • n on F3 (lOe) 

The function fo{x,s) is a given function representing the 
age distribution at initial time t = Q, the function /a is 
used to represent a purely diffusive production of age mass 
along F2 boundary portions, and the total age mass flux is 
J = — 6'DVg. In the Laplace domain the delta distribu- 
tion 5{t) being simply 1, for all times the zero-age boundary 
condition is proportional to the net fluid flux q(a;, t) ■ n. All 
types of boundary conditions can be handled in a standard 
way, just as is done by solute transport codes (given the 
particularities of each of the finite element, finite volume 
or finite differences methods in relation to the implementa- 
tion of boundary conditions). The only difference is that 
boundary functions are expressed with complex numbers. 
However, they still can be treated using the same formalism 
than the one used with real arithmetics functions, because 
in fact the imaginary part of the zero-age condition is always 
zero. The Dirichlet condition is represented by the complex 
function z = 1 + j x Q (where j denotes the imaginary unit 
of complex numbers, = — 1), and the Cauchy condition is 
represented by the complex function z = (\ - n-\- j x Q. Both 
conditions have a zero imaginary part and thus fall into the 
space of real numbers. 

The zero-order source term mo corresponds to an initial 
condition for <; on r for all times t and positions x, and will 
need specific handling since it is not known and not necessar- 
ily zero at any time t, a priori. The only location where mo 
is known is along the inflow boundary Fi, where it is infinity. 
This implies that equation (9) is of nonlinear type because 
it partially depends on the unknown solution g{x, t, 0). The 
term mo should thus not be neglected, in particular at the 
neighborhood of the inflow boundaries, and above all for 
hydrogeological configurations that tend to evolve toward 
increased mixing of ages or that contain important fractions 
of young age masses due to important infiltrations by su- 
perficial recharge. Neglecting this term would exclude the 
correct treatment of the distributions g{x, t, r) for which 



the maximum intensity is located around small values of 
r, as can be found with cases of young and shallow flow 
systems undergoing recharge by rainfall infiltration, or also 
with strongly dispersive systems inducing high mixing of age 
masses. The specific, however still very theoretical case of 
systems described by the exponential age distribution model 
corresponds to the worst configuration for which neglecting 
the term g{x, t, 0) would not allow the solution g to show 
maximum age intensities at exactly r = 0. The numerical 
treatment of this source term will be discussed later in sec- 
tion 3.3. 

Provided an initial distribution g{x, 0, r) and the prescrip- 
tion of appropriate boundary conditions, equation (9) can 
be solved for a finite number of discrete values Sk of the 
Laplace variable s, just as if the system of equations invok- 
ing equation (9) for each discrete value Sk would correspond 
to a multispecies mass transport equation system, with each 
of the species being assimilated to a specific discrete value 
Sk, and within which the species would not interact. By 
letting £fc() = V ■ q() - V ■ 6*0X70 + Sfe6l() -I- go() be the lin- 
ear transport operator associated to variable Sk, the system 
of equations corresponding to Nl discrete Laplace variables 
results in: 

(do \ 
-T^ + Ckj g{x,t,Sk) = rhk , k^l,...,NL (U) 

with rhk — q\gi{sk) + mo. The function gi{sk) is an age 
distribution expressed in the Laplace domain. It is Dirac 
for the classical case (gi — 1) or it needs to be speci- 
fied as gi{x,t,Sk) for cases where internal volumes of wa- 
ter are produced with particular, not necessarily nil age 
distributions. Equation (11) combines a time-evolution of 
Laplace transformed functions, and will consequently be re- 
ferred to as the "Time-Marching Laplace- Transform" ap- 
proach (TMLT). The transport process of complex functions 
operated with the TMLT approach is schematized in Fig- 
ure 2. For fixed x and t, the function ^(r) typically is an 
oscillatory function whose period increases with Sk and de- 
pends on the model coefficients q and D {Sudicky [1989]). 




Figure 2. Schematic representation of the TMLT tech- 
nique used to solve the 5-D age equation. 

The TMLT approach can also be used to solve for con- 
stituent exposure time density only for simplified cases 
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where the advection velocity in the r— direction Va and 
the reaction term require a total independence on the 
independent variable r and on the dependent variable 
Pa{x,t,T). Another acceptable situation is when the prod- 
uct Va{x, T)pa{x, t, t) and the reaction term can both be 
transformed in the Laplace space. Given these restrictions, 
the governing equation to be solved is: 

^+V-(vp„-DVp„) + L{^^,r,s}=r„ (12) 

Possible applications of equation (12) can be met in prob- 
lems defined by an exposure-time velocity Va that fol- 
lows a finite mcmorjr model (for example sec the bacte- 
rial transport case presented by Ginn [2000b]), with which 
Va is represented by a left-continuous Hoavisidc step func- 
tion {H{y) = for J/ < and H{y) = 1 for y > 0), 
Vair) = —H{t), with the consequence that the transfor- 
mation of the term Q-^('Opc«(t) gj^gg _ 

To the author's knowledge, there is no such approach that 
lets time-marching processes operating on complex Laplace 
transformed functions, allowing solution of 5-dimensional 
equations of the form of equation (1). 

3.2. Numerical Inversion of the Transformed Functions 

To return into the real age-domain, the fields g{x, t, s) or 
Pa{x,t,s) need to be numerically inverted. Various authors 
(e.g. see Sudicky [1989]; Cornaton [2003]; Comaton and 
Perrochet [2006a]; Morales- Casique and Neuman [2009]) 
have shown that the numerical inversion of Laplace Trans- 
form finite element solutions arc generally accurate when 
using the algorithm of Crump (Crump [1976]) and when 
further combining it to the quotient-difference algorithm, as 
proposed by de Hoog et al. [1982]. The Laplace Transform 
schemes can suffer from oscillations near discontinuities dur- 
ing the inversion, in relation to the non-uniform convergence 
of the series in the inversion formula. A control of these os- 
cillations can be done by increasing the number of discrete 
Laplace variables Nl — 2n + l and decreasing the mesh/grid 
size near discontinuities. It can be reported that a value for 
n ranging between 10 and 20 generally yields accurate so- 
lutions in the presence of sharp gradients ( Cornaton [2003] ; 
Cornaton and Perrochet [2006a]; Morales- Casique and Neu- 
man [2009]). 

It can also be noted that the TMLT method provides 
not only a solution for the density distribution, but also 
for its integral (i.e. the cumulative distribution G{x,t,T) = 

g{x, t, u)du), when one is taking advantage of the Laplace 
Transform property G{x,t,s) = g{x,t, s) / s. Consequently, 
at any selected time t it is possible to obtain both g{x, s) 
and G{x, s) by applying the inverse transform on g and G. 

3.3. Computational Aspects of the TMLT Methods 

In terms of memory usage, the TMLT method involves 
storage requirements that are proportional to N^NiNl <C 
NxNtNr (Nx, Nt, Nr and A''^ being the number of computa- 
tional points, computational times, age classes and Laplace 
variables, respectively). In most of the cases, the amount 
of data needed to store the age distributions in the complex 
Laplace space will be much smaller than the amount needed 
in the real age domain, rendering the method particularly 
interesting when very old ages are present in the modelled 
system. This amount of discrete data partly depends on the 
shape and complexity of the functions (e.g. multi-modal 
functions would obviously require a larger amount of dis- 
crete time values to allow for a good representation of the 
functions), but also on the maximum age that is to be ac- 
counted for. The bigger the maximum age is, the bigger 
the amount of age classes needed to properly represent the 



age distribution. This limitation is not affecting the TMLT 
method since for any extent of the age dimension the num- 
ber of needed discrete Laplace variables remains unchanged. 

The method is consccjuontly particularly well-suited for real- 
site applications that rc(juirc fine discretization of the com- 
putational mesh/grid. 

The global computational demand of the TMLT method can 
be compared to that of a A^-species mass transport prob- 
lem, with the important simplification that the Nl species 
are independent and do not react. The solving of the A^ 
assembled algebraic systems requires iterative solvers that 
are capable of working with complex-arguments matrices 
and double-precision arithmetics. Even if the memory stor- 
age in the r— axis is highly reduced when working in the 
Laplace domain, the overall computational demand calls for 
parallelism and distributed/shared-memory calculations for 
handling real-site 3-D applications over large time periods. 

The zero-order mass source term mo in equation (9) re- 
quires the knowledge of the spatial distribution of the den- 
sity associated to age r = 0, since this density will not 
be zero for many system configurations, as previously dis- 
cussed in section 3.1. The consequence is that the system of 
Nl equations as formalized with equation (11) is of nonlin- 
ear type. This induces an additional computational burden 
since such a nonlinearity has to be treated using iterative 
techniques such as the Picard or Newton-Raphson methods 
to carry out nonlinear iterations within each time-step At 
to yield solution at time t -\- At. With these classical it- 
erative methods the spatial distribution go = g{T = 0) at 
previous time t is used as initial guess to proceed with a 
series of nonlinear iterations with which the field of go, and 
thus of mo, is updated at each new nonlinear iteration in 
order to obtain the solution at time t + At. The spatial 
distribution of go thus needs to be calculated by applying a 
numerical inverse transform to the function g{s) at the age 
value r = after each nonlinear iteration, and until an ap- 
propriate stopping criterion indicates that convergence has 
been achieved. Convergence criteria will typically corre- 
spond to a threshold limit applied to the residues between 
two nonlinear iterations, and can be formulated e.g. using 
the standard absolute error measure — \ g^^'^^ — Po*'"' ^' li 
where j denotes the iteration level and i the computational 
points index. However, the corresponding computational 
burden remains relatively minor for a priori only the imme- 
diate vicinity of the inflow boundary is concerned. 
As an alternative to the iterative solution, a simplistic ap- 
proach is to assume an explicit scheme for calculating the 
state of the function go at time t + At from its previous state 
at time t, i.e. by letting g{x, t + At, r = 0) ~ g[x, t,T — 0). 
Of course the explicit approach will hardly converge nor 
produce stable results with stiff problems, unless it is asso- 
ciated to small time-steps At to keep the error in the result 
bounded, which in turn may render the solution impractical. 

Nonlinearities such as the ones induced by the depen- 
dence of the flow equation to saturation or fluid density 
are implicitly handled by the TMLT method. For each 
time-level, variable saturation conditions can bo updated, 
as well as fluid density p and viscosity /i, yielding time- and 
space-dependent moisture content 6 = 9{x,t) and velocity 
q = q{9, p, p,t). The age equation problem is then solved 
with the corresponding distribution of 9 and q. Concerning 
the viscosity and density dependency of fluid flow, the cou- 
pled flow and nonlinear thermohaline problems can e.g. be 
solved at a given time-level, preliminary to solving for the 
age problem and go forward to the next time-level. Even 
if in the present work only velocity fields derived from the 
incompressible groundwater flow problem are considered, it 
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can be noted that the method is also suited for operating 
with compressible homogeneous fluid flow problems. In ad- 
dition, we can also argue that the TMLT method does not 
present any restrictive limitation in relation to its adap- 
tation to coupled surface-subsurface flow problems, since 
an integrated and implicit flow solution can simply be per- 
formed at any time-level and the resulting velocity field be 
used to solve for the ago equation at the same time-level. 
Unstable TMLT solutions (e.g. originating from advection- 
domiriatod transport processes) can also bo supported with- 
out any specific restriction by stabilization methods (such as 
Petrov-Galerkin based formulations), or upstream- weighting 
and mass-lumping formulations. 

When discretization techniques are used to solve with 
the TMLT method, the terminology can be refined. For 
instance, when the finite element technique is applied, 
the discretization of the system of equations (11) yields 
the TMLTG technique (" Time-Marching Laplace-Transform 
Galerkin"), in reference to the works of Sudicky [1989] on 
the LTG technique. When the finite- volume technique is ap- 
plied, one may refer to the TMLTFV technique. The finite 
differences technique would equivalently yield the TMLTFD 
technique. 

3.4. The TMLTG Solution Technique 

In this section, we apply the finite element framework 
in order to illustrate a numerical implementation of the 
TMLT technique, and to provide additional details related 
to the computational burden associated to the discretization 
of equation (9). 

3.4.1. Finite Element Formulation 

We consider the specific case of the transient groundwa^ 
ter age distribution equation for one value of the Laplace 
variable Sk'- 



dt 



+ V • {q£ik - OTtVgk) + Ofcfffc = m 



(13) 



with Ofe = go + Sk9 and m = gi + 9g[x, t, 0). Initial condi- 
tions and standard boundary conditions yielding solutions 
of equation (13) arc given by equation (10). 

When a standard Galerkin finite element formulation is 
applied to equation (13), the following discretized equation 
is obtained: 



/ ^C^i^ + ek^k)dQ 

Jn 



dt 



(14) 



/ VN(qV'fc - eD\/i!k)dn 

[ N(6lDVi^fe - qd>k) ■ ndV 

f^iqtpk - eHV-ipk) ■ ndr 



= / NmdQ. 
In 



where the vector N = [Ni,...,Nnn] is the clement shape 
functions vector, that is equal to the weighting functions in 
a standard Galerkin formulation, and where the unknown 
variable is denoted by tpk — ipk{x, t, Sk)- Note that the third 
term in the L.H.S. of equation (14) exhibits the gradient 
projection of the dependent variable along open boundaries 
r+. If not neglected, this term will require specific handling, 
e.g. using the technique proposed by Frind [1988] or Cor- 
naton et al. [2004]. 

On any finite element domain Qe, the unknown variable 
is replaced by a continuous approximation of the form 
^'fe = X]n=i i'k.nNn, whcrc upk.n are the nn nodal unknowns. 
After substitution of this trial solution into (14), and after 



dividing the domain and boundary integrals into picccwise 
elemental contributions over an element Qe, the following 
global matrix system of algebraic equations is obtained by 
summing the ne elemental matrices: 

[Afc]{*,} = {F4-[M]{^} (15) 

where [Ak] is the global (complex) stiffness matrix, [M] the 
mass matrix, and {Fk} the load vector: 

Tie ^ 

[Ak] = (B^DB^ - BqN^ -|- NefeN'^)df2e (16a) 



ne „ 

[m] = E/ 

e=l 



N6IN'' dne 



(16b) 



ne „ 

m = E / 

e=l-'^e 



N(qN' - 6IDB' ) • ndFe + [ NmN' dD.e 

(16c) 

with the gradient matrix B = VN. 

The linear algebraic system (15) needs to be assembled 
and solved Nl times, i.e. for each independent Laplace vari- 
able Sk- However, for one given time- level the mass matrix 
[M] is real and only needs to be assembled once. The com- 
plex stiffness matrix [A^] can be splitted into a real part 
[A''] and a complex part [Ak], with: 



e=l •'^e 



and 



(B6IDB' 



ne p 



BqN^ + NqoN^)df2e 



The real matrix only requires a single evaluation, while the 
complex matrix needs evaluation for each Laplace variable 
Sk- However, the complex elementary matrix [Ak]e is ide- 
ally obtained by scaling the real elementary mass matrix 
[m]e = ISiO'N'^dQe by Sk, [aI]^ = Sfc[M]e. Finally since 
the zero-age boundary condition is always 1 in the Laplace 
space, the load term {F^} remains identical for all the Nj^ 
Laplace variables, and also needs to be assembled only once 
for a given time-level, {Ffc} = {F}. Consequently, the com- 
putational burden for the assembling of the Nl equations 
is relatively low, and once the real-arguments matrices and 
vectors [A''], [M] and {F} have been evaluated for the first 
Laplace variable si, the remaining Nl — 1 assemblings are 
obtained by scaling [M] by Sk and by adding the contribu- 
tion of the accumulation term — [M]{-g^} to {F}. Finally 
note that the matrix operations required for assembling and 
solving the linear algebraic system (16) are the same for 
complex matrices as they are defined for matrices with real 
entries, such that the complexity of (16) does not add any 
particular additional technical difficulty. 
3.4.2. Time Discretization 

The well-known numerical artifacts induced by badly con- 
trolled time-marching schemes affecting standard advection- 
dispcrsion transport solutions also affect the TMLT method 
in a similar way. This implies that the same constraints 
regarding the time-discretization have to be considered to 
avoid numerical instabilities (like the constraints on the 
Peclet and Courant numbers), and that appropriate adap- 
tive time-stepping solution techniques also need to be de- 
signed. The accuracy of the TMLT solutions will depend on 
the adopted time integration scheme and stepsize, exactly 
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as with standard transient advection-dispersion equations. 
Whenever the choice on the stepsize is not appropriate to 
ensure stabihty and/or accuracy, the solutions will have the 

tendency of being affected by numerical diffusion wfien fully 
implicit schemes are used or by numerical oscillations when 
fully explicit or centered schemes are used. 
The time-marching procedure needed to solve the transient 
linear algebraic system (15) can be obtained by a standard 
finite-differences approximation of the time-derivatives. The 
resulting algebraic system can be formalized by the follow- 
ing: 

{*,}*+^* = [C,]-i{*,}* + {QJ* (17) 
with coefficient matrix 

[C,] = (^e[Afc]*+^^' + [M]*+^^*) (18) 
and accumulation term 

{QkY = ((e - l)[Afc]*+^^* + ^[M]*+^^*) {*,}*+{F}*+- 

(19) 

with e varying between and 1. Like with standard trans- 
port scheme, when the explicit scheme is used (e = 0) no 
iterations arc needed since the matrices are evaluated only 
with the initial value {^'j;}*, always constant over a given 
time-step. However, to ensure stability, convergence and ac- 
curacy, the Crank-Nicholson scheme (e = 0.5) is generally 
recommended. The fully-implicit scheme (e = 1) is known to 
yield smoother solutions in relation to its intrinsic diffusive 
nature. Adaptive time-stepping strategies can be designed 
to automatically control the stepsize At. One can classically 
use the evolution of the norm of the complex spectrum g be- 
tween two or more time-steps to predict an optimal stepsize. 



4. Validation and Illustration Examples 

In this section the TMLTG scheme is validated and illus- 
trated by means of analytical and numerical solutions. The 
numerical solutions have been performed with the GW soft- 
ware {Cornaton [2007]), which contains the characteristics of 
the TMLTG method as described in section 3.4. The nonlin- 
ear iterations induced by the zero-order source term mo (a:, t) 
in equation (9) are performed using the Newton-Raphson it- 
erative scheme. Initial conditions for the distribution g arc 
obtained by solving the steady-state equation (8) using the 
LTG formalism {Sudicky [1989]; Cornaton [2003]). 

4.1. One-Dimensional Validation of the TMLTG 
Numerical Solution 

The robustness of the TMLTG scheme is validated us- 
ing the pseudo-analytical 1-D solution derived in the Ap- 
pendix. This analytical solution is taken as reference so- 
lution by considering a finite domain of length L = 200 
meters, discretized into 200 pipe elements of 1 meter length. 



Dispersivity and molecular diffusion are fixed to ai, = 2 m 
and Dm = 0, respectively. An initial equilibrated age dis- 
tribution g{x, 0, r) is considered as being the steady-state 
solution resulting from the uniform and constant velocity 
V = 1 m/d. At t > 0, this velocity is assumed to be two 
times lower, i.e. v = 0.5 m/d. The resulting temporal evo- 
lution of the age distribution g{x, t, r) is shown in Figure 3a 
at the position x — 100 meters. The figure compares the 
analytical and numerical solutions. The TMLTG numerical 
solution makes use of the Crank-Nicholson scheme for the 
time-discretization, with a constant stepsize of 0.05 days. 
The number of Laplace variables used to store the age distri- 
butions is Nl = 2n + 1 = 31. Figure 3b shows the temporal 
evolution of mean age a{x, t) and standard deviation (Jg{x, t) 
at a; = 100 meters. These two quantities have been obtained 
by postprocessing the transient age distribution. Note that 
this is a means for validating the TMLTG solution obtained 
for this example, since the direct mean age solution can be 
taken as a robust reference from a numerical point of view. 
Mean age is two times bigger for the new equilibrium (at 
initial state it equals ^ = ^ = 100 days, and at final state 
it equals ^ = -g^ = 200 days), and so is standard deviation. 
With respect to the evolution in time of the age distribution, 
one can see that from initial time and during approximately 
one pore volume (i.e t = ^ = ^ = 100 days, the turnover 
time corresponding to the initial condition), the initial age 
distribution is simply shifted toward bigger age values, its 
shape remaining intact. After time t = 100 days, the distri- 
bution experiences modifications and displays a successive 
apparition of peaks, a reduction of its maximum intensity 
and the increase of its dispersion. This transition phase is 
occurring until the new equilibrium is met, after about 350 
days (thus more than one initial pore volume plus one new 
pore volume), and it can be described as being the time- 
period during which a crossing of ages occurs. This partic- 
ular time-period corresponds to the phase during which the 
contribution of infiltrated water molecules under the new 
fiow regime to the new state of the age distribution is occur- 
ring. It is clear that this contribution needs, at least, one 
pore volume to be effective. In other words, since velocity 
has been lowered by a factor 2, all water molecules need two 
times the amount of time to travel to a same downstream 
location in comparison with the initial (equilibrated) situa- 
tion. Thus, during at least one (previous) pore volume the 
age distribution is shifted toward older ages. After, it mixes 
with older ages (the ones of molecules that took more than a 
pore volume to reach the same location) and it progressively 
transforms until new equilibrium is met. 
The surprising snail shell shape image displayed in Fig- 
ure 3c corresponds to the time evolution of the complex 
age distribution (still in the Laplace space) at the position 
X = 100 meters. The 2-D graph on the figure plots the imag- 
inary part of g{x, t, s) against its real part, from initial state 
{t = days) to final state (t = 350 days), with a temporal 
resolution of 5 days. The 3-D graph in Figure 3d gives a 
similar space-time representation. 
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Figure 3. Temporal evolution of the 1-D transient-state age distribution g{x, t, r) at a; = 100 m. (a) 
Distribution g{x, t, r) at some selected times; (b) Temporal moments: mean age a{x, t) and standard 
deviation ag{x,t); (c) and (d) Real and imaginary parts of g{s) at some selected times. 



Figure 4 shows the same solutions using 2-D and 3-D 
representations. The shifting phase is well illustrated 
by an evolution of the age distribution in time follow- 
ing an angle of ^ = 22.5 degrees, showing that the in- 
crease of the maximum intensity of ages is of 0.5 per unit 
time, in relation to the new flow regime dictated by a 
two-times smaller velocity as compared to the previous 
flow regime (i.e. there is a ratio of 2 between previous 



and new velocity, or equivalently between the previous 
and new turnover time). This phase is followed by the 
progressive shaping of the new distribution (when ages 
from the previous state are crossed by ages of the new 
state), until the new equilibrium is met. The "crossing 
of ages" phase is characterized by an angle of exactly 
45 degrees, reflecting the fact that all water molecules 
experience the same ageing per unit time. 
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Figure 4. Temporal evolution of the 1-D transient-state age distribution g{x, t, t) at x = 100 m. Left: 
2-D representation; Right: 3-D representation showing the analytical solution in dark grey and the 
numerical solution in grey. 



As discussed in section 3.4.2, the accuracy and phys- 
ical meaning of the TMLTG solutions are directly re- 
lated to the time-marching scheme and the stepsize that 
are employed. We briefly illustrate the effect of the cen- 
tered Crank-Nicholson and implicit schemes on the 1- 
D solution, by making use of some selected stepsizes. 
The results are shown in Figure 5, where the behavior 
of the solutions meets the classical, well-known effects 
of integration schemes on advective-dispersive trans- 
port solutions. With a stepsize of 1 day, the implicit 



scheme generates a very diffusive solution (in relation 
to the numerical diffusion that is artificially added by 
such a scheme) and the Crank-Nicholson scheme gener- 
ates a highly oscillating solution. When the stepsize is 
progressively reduced, the oscillations with the Crank- 
Nicholson scheme are damped and the solution quickly 
converges to the exact solution (met for At = 0.1 days), 
while the implicit scheme becomes less diffusive and is 
getting closer to the exact solution (however not fully 
met for At = 0.01 days). 
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Figure 5. Effects of the time-integration scheme and of the stepsize on the temporal evolution of the 
1-D transient-state age distribution g{x, t, t) at x — 100 m. The figure compares the analytical solution 
to the numerical solutions obtained with the Crank-Nicholson and implicit time-integration schemes, (a) 
Implicit solution; (b) Crank-Nicholson solution. 



4.2. Example of a 2-D Aquitard- Aquifer System with 
Natural Recharge 

In order to illustrate the effect of annual recharge 
fluctuations on age distributions, we make use of a 2-D 
synthetic vertical system that corresponds to an ide- 
alization of a 1000 meters long and 50 meters thick 
aquitard-aquifer system, as shown in Figure 6. The 10- 
meters thick aquitard is located between the elevations 
z = 20 m and z = 30 m. The mesh is discretized using 
bilinear quadrangular elements of size 10 x 1 meters. 
Recharge is supposed to be uniformly distributed on 
top of the model, and is modelled by prescribing fluxes 
(Neumann- type condition). The outlet has a finite size 
of 10 meters, and is located at z = 50 m between the 
points X = 990 m and x = 1000 m. A fixed hydraulic 
head of 51 meters is prescribed. Hydraulic conductiv- 
ity in the two aquifers is K = 10~^ m/s, and aquitard 
hydraulic conductivity is K = 10^^ m/s. The stor- 



age coefficient is S' = 7.5 x 10""^ 1/m in the aquifers 
and S = 10~^ 1/m in the aquitard. Aquifer poros- 
ity is 0.1, and aquitard porosity is 0.35. Dispersivity 
coefficients are nil in the aquitard, while longitudinal 
dispersivity is a/, = 10 m and transverse dispersivity is 
ttT = 1 m in the aquifers. Finally, both aquitard and 
aquifers have the same coefficient of molecular diffusion 
Dm = 2.3 X lO-^m/s^. 

The initial state is obtained by assuming a steady-state 
recharge of intensity 10~^ m/s (the average of the time- 
series applied at transient-state), and is represented 
in Figure 6. The hydraulic solution displays vertical 
drainage effects through the aquitard, downward ori- 
ented between approximately x = and x = 600, and 
upward oriented after x — 600. The steady-state mean 
age distribution is typical of the one occurring in uni- 
formly recharged systems (generating exponential-like 
age distributions), but it also displays the perturbation 
induced by the presence of the aquitard, which gener- 
ates mixing of ages by diffusion. 
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Figure 6. Theoretical 2-D vertical aquitard-aquifer system with uniform recharge, (a) Hydraulic solution 
at steady-state with flow net illustrated by means of the advective tracking of some particles from recharge 
to discharge, (b) Mean age solution at steady state in days. The points PI, P2, P3 and P4 are observation 
points used to show the evolution in time of the hydraulic and age solutions. 



The transient-state scenario corresponds to the ap- 
plication of a cyclic forcing using a recharge rate of 
variable intensity on top of the model, and by keep- 
ing the water level at outlet unchanged. Simulations 
have been performed using the implicit time-integration 
scheme with a constant stepsize of 1 day, and using 
Nl = 2n + 1 = A1 discrete Laplace variables. The 
transient evolution of hydraulic head and mean age at 
the observation points is given in Figures 7a and 7b. 
Mean age variations are high, in relation to the im- 
posed high hydraulic variations. A first observation is 
related to the fact that the steady-state situation (us- 
ing the average of annual fluctuations of infiltration) 
does not provide an initial condition that is close to 
some average, or pseudo-equilibrium, considering that 
transient fluctuations of infiltration are continuously re- 
peated for several years. The new pseudo-equilibrium 
is met only after repeating annual series during at least 
8 years for the aquitard and the lower aquifer, while the 
upper aquifer and the outlet tend to equilibrate around 
annual fluctuations already after 1 year. This observa- 



tion points out the importance of evaluating realistic 
initial conditions for age mass transport problems, the 
steady-state (average) situation being not necessarily 
a proper representation of an initial state. A general 
correlation with the hydrologic events can be done: the 
intuitive tendency of ageing of water during drought 
periods and vice-versa can be observed. This tendency 
however includes some delays between the hydrologic 
and age responses, especially in the isolated, confined 
lower aquifer and in the aquitard. 

The transient age frequency distributions recorded 
at the observation points are given in Figures 7c to 7f 
and in Figure 8. The distributions display a relative 
complex temporal evolution as compared to the smooth 
aspect of the initial, steady-state distributions. All dis- 
tributions become multi-modal with a successive gener- 
ation of several modes (peaks) . The aquitard and lower 
aquifer are invaded by young water due to new infiltra- 
tion rates. This is revealed by the apparition of new 
modes at young ages. As time passes, these new water 
components take age and a displacement of these modes 
towards older age values can be observed, with an angle 
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close to 45 degrees, as a reflection of the unit increase 
of age per unit time (i.e. the line passing by the maxi- 
mum intensity of ages is defined by ^li^i^ ~ l). The 
annual new amount of freshly infiltrated water gener- 
ates new peaks, and the consequence is that age dis- 
tributions become multi-modal. The same observations 
can be made for the points PI and P4 in the upper 
aquifer, with the difference that these points show much 
less tailing effect than points P2 and P3. The linear 
shape of the g-spectrum for points PI and P4 is typi- 
cal of situations that are close to the conditions of the 
perfect exponential mixing (uniform recharge from the 
top along the entire domain). On the opposite, the 



g-spectrum for points P2 and P3 differs from the lin- 
ear, exponential-like case, and is more similar to that 
of the 1-D example presented in section 4.1 (see Fig- 
ure 3). An important observation is that, as opposed 
to the steady-state situation for points PI and P4 (that 
both display the theoretical situation of an exponen- 
tial mixing), the transient response does not present 
the same behavior but rather displays age distributions 
characterized by multi-modality and values at the ori- 
gin that strictly differ from the exponential case (with 
which gij — Q) ~ I/tq, with tq denoting the turnover 
time, or ratio between porous volume and steady-state 
discharge rate). 




Figure 7. Temporal evolution of flow, mean age and age distribution in the theoretical 2-D vertical 
aquitard-aquifer system of Figure 6. (a) Hydraulic head; (b) Mean age; (c) Age distribution above the 
aquitard; (d) Age distribution at the outlet; (e) Age distribution within the aquitard; (f) Age distribution 
under the aquitard. 
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Figure 8. Temporal evolution of the age frequency distributions in the theoretical 2-D vertical aquitard- 
aquifer system of Figure 6, using a surface representation. 



4.3. Example of a Pumping- Well in a 3-D Homogeneous 
Aquifer with Natural Recharge 

Our last example corresponds to the situation of 
the activation of a pumping- well inside a homogeneous 
aquiferous block of 100 by 100 meters and of thickness 
50 meters. The initial state is obtained by recharg- 
ing the aquifer at steady-state at a uniform infiltration 
rate 2 x 10^* m/s on top of the model. The outlet 
is modelled using a fixed hydraulic head of 50 meters 
prescribed along one vertical face, as indicated in Fig- 
ure 9a. As opposed to the previous 2-D example, we 
now make use of a natural, undisturbed initial condi- 
tion for the transient solution. The mesh is made of 
8-noded brick elements of size 1x1x1 meter. Hy- 
draulic conductivity is set to K — 10^^ m/s, storage 
coefficient is S* = 3 x 10"^ 1/m, porosity is 0.1, longi- 
tudinal dispersivity is = 1 ni, transverse horizontal 
dispersivity is ar./i = 0.1 m, transverse vertical disper- 



sivity is aT,v = 0.01 m, and the coefficient of molecular 
diffusion is — 2.3 x lO^^m/s^. 
The transient-state situation starts with the activation 
of a pumping-well located at {x,y,z} — {50,50,25}, 
using an extraction function that globally corresponds 
to cycles of pumping during 1 week at a rate rang- 
ing between 2 x lO^^m'^/s and 1.5 x 10~^m^/s, and 
no-pumping during 1 week. The remaining boundary 
conditions are kept unchanged. Figure 9a shows the 
hydraulic solution at a given simulation time, using 
three isosurfaces of hydraulic head. Figure 9b shows 
the mean age solution. Mean age tends to meet pseudo- 
equilibrium (i.e. an average situation considering the 
pump rate fluctuations) after 3 to 4 years. The solu- 
tion displays important variations in the neighborhood 
of the well. At the well location, mean age slightly 
increases due to the capture of water molecules that 
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have traveled from upstream locations and that enrich 
the mixing process with older ages. The small increase 
(from 5.75 years at initial time to around 6 years) at 
the well may be explained by the result of a mixing 
that equilibrates a component of newly captured older 
ages and a component of newly captured young ages 
coming from top locations above the well. The tempo- 
ral variations of mean age due to pump rate variations 
are of about 0.5 year. Only 2 meters above the well lo- 
cation, mean age strongly decreases (from 5.45 years at 
initial time to around 2.5 years, i.e. more than 50% or 
relative change) due to the capture of a major compo- 
nent of young water that infiltrates above the well. The 
temporal variations of mean age are of about 0.75 year. 
These important changes in mean age point out the 
risk associated to the interpretation of apparent ages 



Initial di! 
{no PL 



5. Conclusions and Outlooks 

A novel and promising numerical scheme for solving 
the water age distribution problem in environmental 



deduced from isotopic methods, particularly when in- 
terpretation is carried out for estimating flow velocities 
or simply constraining the calibration of a flow model. 
The age distribution solution makes use of = 2n + 
1 = 41 discrete Laplace variables and of the implicit 
time-integration scheme with a constant stepsize of 1 
day. The effect of activating the pump on the age dis- 
tribution recorded at the well location is shown in Fig- 
ure 9c. The shape of the distributions confirms the 
interpretation based on the analysis of mean age solu- 
tion. The distribution at the well node is more dispersed 
than the distribution recorded 2 meters above, and dis- 
plays a bigger contribution of young water. The pump 
fluctuations are well-marked and revealed by peak vari- 
ations. 



transient flow domains of arbitrary complexity has been 
elaborated. The presented scheme can be of high in- 
terest for various interconnected fields like atmospheric 
and ocean circulation modelling or surface and sub- 
surface hydrology. The main characteristic of the al- 
gorithm is that a reduction in problem dimension is 



(a) . (b) 




25 



Figure 9. Temporal evolution of mean age and of the age frequency distribution at the pumping-welL 
(a) Hydraulic solution at a given simulation time illustrating the effect of the pumping; (b) Mean age 
evolution; (c) Age distribution evolution; (d) Age distribution evolution in the complex Laplace space. 
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obtained by mathematical transformation of the age 

dimension using the Laplaee transform operator. A 
standard time-marching procedure is then apphed to 
the Laplace transformed age distributions to simulate 
their temporal evolution. The main advantage of trans- 
forming the age distribution function is to be relieved 
from the need of discretization in the age-axis. The 
age dimension is therefore compressed in the complex 
Laplace space, which is independent on the age mag- 
nitude. A direct consequence is that, as opposed to 
the few existing methods, age distributions that dis- 
play very big age values (e.g. of the order of 1 Million 
years or more, as observed in deep sedimentary basins) 
can be modelled without any specific limitation. The 
proposed algorithm is consequently particularly well- 
suited to treat hydrologic or atmospheric problems in- 
volving predictions of circulations and mass transfers 
over very large time-scales (in both the clock-time and 
age coordinates). Moreover, the algorithm allows for 
the integration of nonlinearities affecting a flow solu- 
tion, such as the ones induced by variably saturated 
conditions in subsurface hydrology, surface-subsurface 
coupling of flow dynamics, or buoyancy effects encoun- 
tered in thermohaline processes, and can also be ap- 
plied to compressible homogeneous fluid flow problems. 
The approach can be implemented using the classical 
numerical integration schemes like the finite element, 
finite volume or finite differences methods. 
Since the method combines a Laplace transform of 
the age equation in the age-coordinate to a classi- 
cal finite-difference approximation of the time evolu- 
tion in the clock-time coordinate, the well-known nu- 
merical artifacts generated by badly controlled time- 
marching schemes and affecting the standard advection- 
dispersion transport equations also need to be consid- 
ered. Consequently, the constraints on the space- and 
time-discretization have to be considered to avoid nu- 
merical instabilities (like the constraints on the Peclet 
and Courant numbers). 

The approach computational demand can be compared 
to that required by the solution of a transport prob- 
lem involving independent species with first-order de- 
cay (however without reactions between the species), 
the number of species being equal to the number of 
used Laplace variables (ideally between 10 and 20, at 
maximum 40). Thus in order to be of real practical 
interest for handling problems of big size (i.e. with 
meshes/grids of more than 1 to 10 Million unknowns) 
and high parameter contrasts, the need of paralleliza- 
tion techniques applied to the assembling and solving 
parts of the solution can be foreseen. 
The presented numerical experiments do not have the 
pretention of providing an exhaustive analysis of the 
effects of time-varying flow conditions on age distribu- 
tions. They are rather used to validate and illustrate 
the proposed method. They however show that tran- 
sient age distributions undergoing variations in recharge 
and/or discharge conditions can reveal to be very com- 



plex and that realistic initial conditions for age can be a 

critical aspect to be resolved. These experiments point 
out probable severe complications in the field of age dat- 
ing with single and multiple-tracer techniques and flow 
model calibration based on isotopic interpretations of 
age. Theoretical analysis of the time-dependent behav- 
ior of age distributions should be carried out in order to 
gain insights on fluid and mass transfer processes over 
various time-scales and space-scales. 
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Appendix: Derivation of a 1-D Pseudo- 
Analytical Solution for Transient Age 

We consider a semi-infinite domain (0 < a; < oo) witli 
uniform fluid velocity v. At any position x, we suppose 



tliat tlie initial age distribution is of Gaussian type, e.g. 
as a rcstilt of an equilibrated distribution in a uniform 
and constant velocity field u: 

X / {x — utY" 



9u{x,t) 



(Al) 



where D = ulU -\- Dm- Tlie transient-state age distri- 
bution g ~ g{x,t,T) is solution of the 1-D version of 
equation (4): 



dg^ ^ _^^dg_ ^ ^cfg^ _ dg^ 
dt dx dx"^ Dt 



(A2) 



in which the dispersion coefficient reads D = aLV-\-Dm, 
and with the initial and boundary conditions: 



g{x,0,T) = gu{x,T) 

g{0,t,T) = d{T) 
,dg 



- D- 



dx 



= 



(A3) 
(A4) 
(A5) 



To obtain a solution to that problem, we first apply 
a Laplace transform to equations (A2) to (A5) in the 
T-dimension: 



9g 
dt 



dg d g 



dx 



D 



dx'^ 



g{x,0,s) =gu{x,s) = exp 



- D 



g{0,t,s) = l 
dg 



sg-\-g{x,t,0) (A6a) 
iDsW 



(A6b) 
(A6c) 




dx 



= 



(A6d) 



where s is the complex Laplace variable dual to age 
T, and where g — g{x,t,s) is the transformed state of 
g{x, t, t). Now applying a second Laplace transform to 
equations (A6a), (A6c) and (A6d) in the f-dimension 
yields: 



— V 



dg , ^d^g 



dx 



D 



dx"^ 



{r-\- s)g + gu = Q 
1 



g(0,r,s) = 



- D 



di 
dx 



= 



(A7a) 
(A7b) 
(A7c) 



where r is the complex Laplace variable dual to t, 

g = g{x, r, s) is the transformed state of g(x, t, s), and 
where we assumed that g(x,t,Q) = 0. This assumption 
is valid since for this specific configuration of a flow sys- 
tem, the frequency for having zero ages is not nil only at 
a; = where it is Dirac. The solution to this boundary 
value problem reads: 



x(l-,8„(^+«)) 



g{x,r,s) = 



re 



■ sue 



r{r -\- su)) 



with 



w = 1 - - 
u 



(A8) 



(A9a) 
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iaLZ 



AaLZ 



(A9b) 



(A9c) 



and whore molecular diffusion has been neglected 
(-Dm = 0) for a sake of simplicity. The coefficient w 
represents the relative change in velocity. Inversion of 
equation (A8), g{x,t,s) = L~^{g{x,r,s),r,t}, yields: 

g{x,t,s) = %{x,t,s) + 'ji{x, t,s) - ^2{x,t,s) (AlO) 

with 



lo{x,t,s) = e" 

1 x(l-g„(a)) 



crfc 



1 x(l+/3„(s)) 

+ -e erfc 



^2{x,t,s) = -■yo{x,t,s)erfc 



X — Py{s)vt 

x + I3y{s)vt 
/ X - /3u{s)vt 



(Alia) 

(Allb) 



(Allc) 



1 a!(l+g„(s)) 



■suit 



erfc 



x + l3u{s)vt 



2^/aLvt 

One can verify that the mass in the system is al- 
ways 1 by evaluating the zero-order age moment of 
g{x,t,T), g{x,t,T)dT = g{x,t,s ^ 0) = 1. The 
solution (AlO) can be expressed in dimensionless form 
by letting L be a characteristic length, and by us- 
ing the change of variables X = j^, T = ^ = ^, 
Pe = ^ = — , and z = TyS. Note that Ty = - denotes 
the mean turnover time (or reservoir mean residence 
time) of tho new equilibrated flow system (under the 
flow regime associated to velocity v) , and that Pe is the 
Peclet number. The dimensionless solution becomes: 

g{X, T, z) = 70 (^, T, z) + 71 (X, T, z) - 72 (X, T, z) 

(A12) 

with 



MX,T,z) 



XPe(l-Dy,(z)) 



■zojT 



(A13a) 



1 fx 

72(X,T,z) = -7o(^,T,z)erfc| 



1 XPe(l + fitu(z)) 



\^ 2,/TjP^ J 

(A13c) 

X + /3Uz)t" 



-zujT 



erfc 



2y/T/P~e 



and we have g{x, t, s) = Tyg{X, T, z). When time tends 
toward infinity, the new equilibrated age distribution 
9v{X, C) = g{X, 00, C) = L~^{g{X, 00, z), z, Q is found: 



9v{X,0 



XVPe 



47^(3/2 



exp 



4C 



(AM) 



with the dimensionless age dimension variable ( = ^. 
The first-order age moment of equation (AlO) gives the 
transient-state mean age a{x,t): 



a{x,t) 



lim 

s=0 



dg{x,t,s) 



ds 



(A15) 
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- +wt 
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+ —{x — vt)evk 
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X — vt 



/ x + vt\ 



which satisfies the properties a (a;, 0) = ^ and a{x,(x) = 
^ , and that is solution of the 1-D form of the mean age 



equation (see equation (7a)), || 



dx 



when the boundary condition a(0,t) = is used. Di- 
mensionless mean age finally reads: 



a{X,T) = {l-w)X + wT 



(A16) 



+ -(X-T)erfc 
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